“Scrap_Proyecto_TS”
Alguna idea de cómo podemos ir trabajando y puliendo el proyecto.
# Gestion de Datos
#Leemos archivo
datos.R <- read_csv('./Series_Remesas_Proyecto_Team3.csv')## Parsed with column specification:
## cols(
## Fecha = col_character(),
## `Remesas Familiares Total Millones USD` = col_double(),
## `Remesas Familiares Money Orders Millones USD` = col_double(),
## `Remesas Familiares Cheques Personales Millones USD` = col_double(),
## `Remesas Familiares Transferencias Electronicas Millones USD` = col_double(),
## `Remesas Familiares Efectivo y Especie Millones USD` = col_double(),
## `Remesas Familiares Total Miles de Operaciones` = col_double(),
## `Remesas Familiares Money Orders Miles de Operaciones` = col_double(),
## `Remesas Familiares Cheques Personales Miles de Operaciones` = col_double(),
## `Remesas Familiares Transferencias Electronicas Miles de Operaciones` = col_double(),
## `Remesas Familiares Efectivo y Especie Miles de Operaciones` = col_double()
## )datos.R$Fecha <- as.Date(datos.R$Fecha, format = '%d / %m / %y')
###########################################
# Selección de posibles modelos
# Graficamos la serie de tiempo sin modificaciones
TS <- ts(datos.R[, 2], start = 1996, frequency = 12)
TS <- as.quarterly(TS, FUN = sum)
autoplot(TS)
FAC <- acf2(TS, 24)
adf.test(TS)## Warning in adf.test(TS): p-value greater than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: TS
## Dickey-Fuller = -0.18868, Lag order = 4, p-value = 0.99
## alternative hypothesis: stationarybestNormalize::boxcox(TS)## Standardized Box Cox Transformation with 103 nonmissing obs.:
## Estimated statistics:
## - lambda = 0.6377471
## - mean (before standardization) = 354.7057
## - sd (before standardization) = 131.704# Mientras que la estimación de lamda para la transformación Box-Cox sugiere una
#transformar la serie con una transformación Box-Cox con lambda = 0.63,
#favoreciendo la interpretabilidad del modelo y considerando que el cambio de
#varianza no es tan extremo se tomará la transformación log()
# Notamos algo de incremento de varianza mientras aumenta el nivel,
#por lo tanto aplicamos log para estabilizar la varianza
log_TS <- log(TS)
autoplot(log_TS)
FAC <- acf2(log_TS, 24)
adf.test(log_TS)##
## Augmented Dickey-Fuller Test
##
## data: log_TS
## Dickey-Fuller = -1.8394, Lag order = 4, p-value = 0.6433
## alternative hypothesis: stationary# La serie claramente no es estacionaria, muestra un aumento constante de nivel
#y un valor-p de 0.64 en la prueba de Dickey-Fuller por lo tanto aplicamos
#una diferencia.
D1_log_TS <- diff(log_TS)
autoplot(D1_log_TS)
FAC <- acf2(D1_log_TS, 24)
adf.test(D1_log_TS)##
## Augmented Dickey-Fuller Test
##
## data: D1_log_TS
## Dickey-Fuller = -2.6152, Lag order = 4, p-value = 0.3216
## alternative hypothesis: stationary# La serie parecen tener estacionalidad anual ya que todos los múltiplos de 4
#en la FAC tienen un pico, reemplazamos nuestra diferencia normal por una
#estacional con un lag de 4.
ggseasonplot(D1_log_TS)
SD4_log_TS <- diff(log_TS, 4)
autoplot(SD4_log_TS)
FAC <- acf2(SD4_log_TS, 24)
adf.test(SD4_log_TS)##
## Augmented Dickey-Fuller Test
##
## data: SD4_log_TS
## Dickey-Fuller = -1.9137, Lag order = 4, p-value = 0.6123
## alternative hypothesis: stationary# Mientras que se ha tomado en cuenta el efecto de la estacionalidad,
#una prueba Dickey-Fuller con valor-p de 0.61 y la FAC nos dejan claro que la
#serie aún no es estacionaria. Por lo tanto aplicamos una diferencia no
#estacional.
D1_SD4_log_TS <- diff(SD4_log_TS, 4)
autoplot(D1_SD4_log_TS)
FAC <- acf2(D1_SD4_log_TS, 24)
adf.test(D1_SD4_log_TS)## Warning in adf.test(D1_SD4_log_TS): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: D1_SD4_log_TS
## Dickey-Fuller = -4.2702, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary# Vemos que todavía hay algunos comportamientos que eliminar de la FAC para
#que parezca ruido.
# El primero que buscamos eliminar es el pico negativo en el cuarto valor
#tanto de la FAC como de la FACP, para esto podemos agregar un término AR
#estacional o un término MA estacional. Dado que este pico sólo se ve en el
#primer múltiplo de 4 en la FAC y luego se ve pero a mucho menor escala, es
#más probable que el término sea MA.
fit1 <- Arima(
log_TS,
order = c(0, 1, 0),
seasonal = list(order = c(1, 1, 0), period = 4),
lambda = NULL,
include.constant = TRUE
)
mean(fit1$residuals)## [1] 0.001016048# SUPUESTO 1: Tienen media cero los residuos
checkresiduals(fit1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,0)(1,1,0)[4]
## Q* = 18.661, df = 7, p-value = 0.00932
##
## Model df: 1. Total lags used: 8# SUPUESTO 2: Tiene un ligero decremento en varianza, nada demasiado grande
Box.test(fit1$residuals,
type = "Ljung",
lag = 24,
fitdf = 0)##
## Box-Ljung test
##
## data: fit1$residuals
## X-squared = 35.826, df = 24, p-value = 0.05706# SUPUESTO 3: Hay suficiente evidencia para decir que los errores están
#correlacionados, se rompe el modelo (con un 90% de certeza)
3 * sqrt(var(fit1$residuals))## [1] 0.2199538# SUPUESTO 5: Parece haber 2 observaciones aberrantes fuera de
#+-3*desviación estándar
coeftest(fit1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## sar1 -0.446877 0.089466 -4.9949 5.885e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# SUPUESTO 6: No hay un modelo con menos parámetros y por lo tanto es
#parsimonioso, el parámetro es muy significativo
autoplot(fit1)
# SUPUESTO 7: El modelo es admisible, sus raíces están dentro del
#círculo unitario
# SUPUESTO 8: Dado que hay sólo un parámetro, no hay correlación entre
#parámetros.
FAC <- acf2(fit1$residuals)
summary(fit1)## Series: log_TS
## ARIMA(0,1,0)(1,1,0)[4]
##
## Coefficients:
## sar1
## -0.4469
## s.e. 0.0895
##
## sigma^2 estimated as 0.005654: log likelihood=114.62
## AIC=-225.24 AICc=-225.11 BIC=-220.07
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.001016048 0.07296822 0.05304835 0.01048119 0.634623 0.4299193
## ACF1
## Training set -0.2290974# AIC = -225.24
# Veamos ahora un parámetro MA estacional
fit2 <- Arima(
log_TS,
order = c(0, 1, 0),
seasonal = list(order = c(0, 1, 1), period = 4),
lambda = NULL,
include.constant = TRUE
)
mean(fit2$residuals)## [1] -0.002683247# SUPUESTO 1: Tienen media cero los residuos
checkresiduals(fit2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,0)(0,1,1)[4]
## Q* = 8.8841, df = 7, p-value = 0.2611
##
## Model df: 1. Total lags used: 8# SUPUESTO 2: Tiene un ligero decremento en varianza, nada demasiado grande
Box.test(fit2$residuals,
type = "Ljung",
lag = 24,
fitdf = 0)##
## Box-Ljung test
##
## data: fit2$residuals
## X-squared = 26.442, df = 24, p-value = 0.3311# SUPUESTO 3: No hay suficiente evidencia para decir que hay autocorrelación
#importante los residuos
# SUPUESTO 4: Se asemejan a una distribución normal los residuos
3 * sqrt(var(fit2$residuals))## [1] 0.1962862# SUPUESTO 5: Parece haber 2 observaciones aberrantes fuera de +-3*desviación
#estándar, las mismas que en el modelo anterior.
coeftest(fit2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## sma1 -0.881435 0.097334 -9.0557 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# SUPUESTO 6: No hay un modelo con menos parámetros y por lo tanto es
#parsimonioso, el parámetro es muy significativo
autoplot(fit2)
# SUPUESTO 7: El modelo es admisible, sus raíces están dentro del círculo
#unitario
# SUPUESTO 8: Dado que hay sólo un parámetro, no hay correlación entre parámetros.
FAC <- acf2(fit2$residuals)
summary(fit2)## Series: log_TS
## ARIMA(0,1,0)(0,1,1)[4]
##
## Coefficients:
## sma1
## -0.8814
## s.e. 0.0973
##
## sigma^2 estimated as 0.004509: log likelihood=123.16
## AIC=-242.31 AICc=-242.18 BIC=-237.14
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.002683247 0.06516561 0.04804542 -0.03127341 0.5751962 0.3893741
## ACF1
## Training set -0.1087152# AIC = -242.31
# La FAC de los residuos del primer modelo indica que podría beneficiarse de
#un término MA no estacional.
fit3 <- Arima(
log_TS,
order = c(0, 1, 1),
seasonal = list(order = c(1, 1, 0), period = 4),
lambda = NULL,
include.constant = TRUE
)
mean(fit3$residuals)## [1] 0.00123736# SUPUESTO 1: Tienen media cero los residuos
checkresiduals(fit3)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(1,1,0)[4]
## Q* = 9.798, df = 6, p-value = 0.1334
##
## Model df: 2. Total lags used: 8# SUPUESTO 2: Tiene un ligero decremento en varianza, nada demasiado grande
Box.test(fit3$residuals,
type = "Ljung",
lag = 24,
fitdf = 0)##
## Box-Ljung test
##
## data: fit3$residuals
## X-squared = 26.223, df = 24, p-value = 0.342# SUPUESTO 3: No hay suficiente evidencia para decir que hay autocorrelación
#importante los residuos
# SUPUESTO 4: Se asemejan a una distribución normal los residuos
3 * sqrt(var(fit3$residuals))## [1] 0.2144822# SUPUESTO 5: Hay sólo una observación fuera de +-3*desviación estándar
#de los residuos y otra muy cerca de estarlo
coeftest(fit3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.218837 0.098777 -2.2155 0.02673 *
## sar1 -0.448765 0.089412 -5.0191 5.192e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# SUPUESTO 6: El parámetro MA estacional parece no ser tan significativo
#como los parámetros de los otros dos modelos, aunque sigue siendo bastante
#significativo. El parámetro AR estacional sigue siendo muy significativo.
autoplot(fit3)
# SUPUESTO 7: El modelo es admisible, sus raíces están dentro del círculo
#unitario
FAC <- acf2(fit3$residuals)
summary(fit3)## Series: log_TS
## ARIMA(0,1,1)(1,1,0)[4]
##
## Coefficients:
## ma1 sar1
## -0.2188 -0.4488
## s.e. 0.0988 0.0894
##
## sigma^2 estimated as 0.005433: log likelihood=117.05
## AIC=-228.11 AICc=-227.85 BIC=-220.35
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.00123736 0.07115693 0.05263366 0.01323543 0.6303178 0.4265585
## ACF1
## Training set -0.009098219# AIC = -228.11El primer modelo viola el supuesto de autocorrelación de los errores y por lo tanto lo desechamos por completo. Nos queda decidir entre el segundo y tercer modelo. La validación de los supuestos entre el segundo y tercer modelo son muy similares. Hay el segundo modelo tiene 2 observaciones aberrantes mientras que el tercero sólo una, pero está muy cerca de tener una segunda y este no es un factor tan importante del modelo. Podemos decir que el segundo modelo es más parsimonioso que el tercero ya que el agregar un parámetro no estacional no tan significativo se complica un poco el modelo. Este tampoco es un punto tan importante ya que aunque no tanto como el parámetro estacional, sigue siento significativo el parámetro MA no estacional agregado. Por lo tanto nos queda el AIC que es mejor en el segundo modelo que en el tercero, esta siendo la diferencia más significativa entre ambos y por la cual escogeremos el segundo.
ARIMA_max <- auto.arima(log_TS, trace = TRUE)##
## ARIMA(2,1,2)(1,1,1)[4] : Inf
## ARIMA(0,1,0)(0,1,0)[4] : -205.428
## ARIMA(1,1,0)(1,1,0)[4] : -228.1737
## ARIMA(0,1,1)(0,1,1)[4] : -241.03
## ARIMA(0,1,1)(0,1,0)[4] : -208.0298
## ARIMA(0,1,1)(1,1,1)[4] : Inf
## ARIMA(0,1,1)(0,1,2)[4] : -241.154
## ARIMA(0,1,1)(1,1,2)[4] : Inf
## ARIMA(0,1,0)(0,1,2)[4] : -242.1223
## ARIMA(0,1,0)(0,1,1)[4] : -242.1839
## ARIMA(0,1,0)(1,1,1)[4] : Inf
## ARIMA(0,1,0)(1,1,0)[4] : -225.1108
## ARIMA(0,1,0)(1,1,2)[4] : Inf
## ARIMA(1,1,0)(0,1,1)[4] : -241.1621
## ARIMA(1,1,1)(0,1,1)[4] : -239.3051
##
## Best model: ARIMA(0,1,0)(0,1,1)[4]ARIMA_max## Series: log_TS
## ARIMA(0,1,0)(0,1,1)[4]
##
## Coefficients:
## sma1
## -0.8814
## s.e. 0.0973
##
## sigma^2 estimated as 0.004509: log likelihood=123.16
## AIC=-242.31 AICc=-242.18 BIC=-237.14# La función auto.arima valida nuestra elección de modelo, considerando al
#SARIMA(0,1,0)(0,1,1)[4] la mejor elección.
# Veamos el modelo graficado contra el log de la serie
autoplot(log_TS) + autolayer(fit2$fitted, series = "SARIMA(0,1,0)(0,1,1)[4]")
# Y graficado con la serie original
autoplot(TS) + autolayer(exp(fit2$fitted), series = "SARIMA(0,1,0)(0,1,1)[4]")
# Veamos el forecast 3 años adelante para el log de la serie
pred_log_TS <- forecast(fit2, h = 12)
plot(pred_log_TS)
# Veamos el forecast 3 años adelante para la serie original.
pred_TS = forecast(fit2, h = 12)
pred_TS$mean = exp(pred_TS$mean)
pred_TS$lower = exp(pred_TS$lower)
pred_TS$upper = exp(pred_TS$upper)
pred_TS$x = exp(pred_TS$x)
pred_TS$fitted = exp(pred_TS$fitted)
pred_TS$residuals = exp(pred_TS$residuals)
plot(pred_TS)
#sarima_forecast <- sarima.for(log_TS, n.ahead=12,p=0,d=1,q=1,P=1,D=1,Q=0,S=4)